home *** CD-ROM | disk | FTP | other *** search
/ TeX 1995 July / TeX CD-ROM July 1995 (Disc 1)(Walnut Creek)(1995).ISO / graphics / mfpic / alfatest / graphbase.mf < prev    next >
Text File  |  1994-07-20  |  19KB  |  1,053 lines
  1. %%%
  2. %%%  File: graphbase.mf
  3. %%%
  4.  
  5. mode_setup;
  6. message "mfpic version 0.2.9 alpha  Tue 21 July 1994";
  7.  
  8.  
  9. %% Global Variables.
  10.  
  11. % when the numeric variable `graphbase' is known,
  12. % then graphbase.mf has been input.
  13. % (idea taken from DEK's cmbase.mf)
  14.  
  15. if known graphbase:
  16.   errmessage "You have loaded graphbase.mf more than once!";
  17. else:
  18.   numeric graphbase;
  19.   graphbase := 1;
  20. fi
  21.  
  22. % pen width
  23. % in pixel coordinates
  24.  
  25. newinternal penwd;
  26. interim penwd:=0.5pt;
  27.  
  28. % temporary path
  29.  
  30. path trash;
  31.  
  32.  
  33. %% Utility Macros.
  34.  
  35. % convert text pairs list t to array p.
  36. % best used inside a vardef.
  37.  
  38. def textpairs (suffix p) (text t) =
  39.  save p;
  40.  pair p[];
  41.  p:=0;
  42.  for a=t:
  43.   p[incr p]:=a;
  44.  endfor
  45. enddef;
  46.  
  47. vardef floorpair(expr p) =
  48.  (floor (xpart p), floor (ypart p))
  49. enddef;
  50.  
  51. vardef ceilingpair(expr p) =
  52.  (ceiling (xpart p), ceiling (ypart p))
  53. enddef;
  54.  
  55. vardef hroundpair(expr p) =
  56.  (hround (xpart p), hround (ypart p))
  57. enddef;
  58.  
  59. vardef minpair(expr u)(text t) =
  60.  pair p;  numeric x, y;
  61.  p:=u;
  62.  for q=t:
  63.   x:=min(xpart p, xpart q);
  64.   y:=min(ypart p, ypart q);
  65.   p:=(x,y);
  66.  endfor;
  67.  p
  68. enddef;
  69.  
  70. vardef maxpair(expr u)(text t) =
  71.  pair p;  numeric x, y;
  72.  p:=u;
  73.  for q=t:
  74.   x:=max(xpart p, xpart q);
  75.   y:=max(ypart p, ypart q);
  76.   p:=(x,y);
  77.  endfor;
  78.  p
  79. enddef;
  80.  
  81.  
  82. %% Coordinate Conversion.
  83. % (affine transformation)
  84. % graph -> pixel
  85. % beginchar defines w and h
  86.  
  87. def xconv(expr xvalue) =
  88.  ((xvalue-xneg)/(xpos-xneg))*w
  89. enddef;
  90.  
  91. def yconv(expr yvalue) =
  92.  ((yvalue-yneg)/(ypos-yneg))*h
  93. enddef;
  94.  
  95. transform ztr;
  96.  
  97. def setztr =
  98.  ztr:=identity
  99.  shifted -(xneg,yneg)
  100.  xscaled (w/(xpos-xneg))
  101.  yscaled (h/(ypos-yneg));
  102. enddef;
  103.  
  104. def zconv (expr a) =
  105.  a transformed ztr
  106. enddef;
  107.  
  108.  
  109. %% Initial Setup.
  110.  
  111. % active_plane is the active drawing plane.
  112. % currentpicture is unknown at this stage,
  113. % so use a def, not a picture assignment.
  114.  
  115. def active_plane = currentpicture enddef;
  116.  
  117. def initpic =
  118.  setztr;
  119.  def active_plane = currentpicture enddef;
  120.  pickup pencircle scaled penwd;
  121. enddef;
  122.  
  123.  
  124. %% Extra Trigonometric and Hyperbolic Functions.
  125.  
  126. def acos (expr x) = angle ((x, 1+-+x)) enddef;
  127. def asin (expr y) = angle ((1+-+y, y)) enddef;
  128.  
  129. def exp primary X = (mexp (256 * X)) enddef;
  130. def ln  primary X = (mlog (X) / 256) enddef;
  131.  
  132. vardef cosh primary X =
  133.   save temp;
  134.   temp = exp X;
  135.   (temp + 1/temp) / 2;
  136. enddef;
  137.  
  138. vardef sinh primary X =
  139.   save temp;
  140.   temp = exp X;
  141.   (temp - 1/temp) / 2;
  142. enddef;
  143.  
  144. def acosh (expr y) = ln (y + (y+-+1)) enddef;
  145. def asinh (expr y) = ln (y + (y++1)) enddef;
  146.  
  147.  
  148. %% Coordinate Systems and Transformations.
  149.  
  150. % Coordinate nesting.
  151.  
  152. def bcoords =
  153.   begingroup
  154.     save old_t;
  155.     transform old_t;
  156.     old_t = currenttransform;
  157. enddef;
  158.  
  159. def ecoords =
  160.     currenttransform := old_t;
  161.   endgroup
  162. enddef;
  163.  
  164. % Coordinate changes.
  165. % Example:  `apply_t (rotated) (theta);'.
  166.  
  167. def apply_t (text transformer) =
  168.   currenttransform := identity transformer
  169.     transformed currenttransform;
  170. enddef;
  171.  
  172. let xslanted = slanted;  % (x+sy, y).
  173.  
  174. def yslanted primary s =  % (x, y+sx).
  175.   transformed
  176.    begingroup
  177.      transform T_;
  178.      origin transformed T_ = origin;
  179.      (1,0)  transformed T_ = (1,s);
  180.      (0,1)  transformed T_ = (0,1);
  181.      T_
  182.    endgroup
  183. enddef;
  184.  
  185. def zslanted primary p =  % (xu+yv, xv+yu), where p = (u,v).
  186.   transformed
  187.     begingroup
  188.       transform T_;
  189.       xpart T_ = ypart T_ = 0;
  190.       xxpart T_ = yypart T_ = xpart p;
  191.       xypart T_ = yxpart T_ = ypart p;
  192.       T_
  193.     endgroup
  194. enddef;
  195.  
  196. def xyswapped =
  197.   zslanted (0,1)
  198. enddef;
  199.  
  200. def boosted primary X =  % boosts for special relativity
  201.   zslanted (cosh X, sinh X)
  202. enddef;
  203.  
  204. % Rotate path f about point p by angle th,
  205. % where f is in _pixel_ coordinates,
  206. % and p and th are in _graph_ coordinates.
  207.  
  208. vardef rotatedpath(expr p,th) expr f =
  209.  f transformed inverse ztr rotatedaround (p,th) transformed ztr
  210. enddef;
  211.  
  212.  
  213. %% Bitmaps, Clipping and Rendering.
  214.  
  215. % bitmap to bitmap --- bitwise operations.
  216.  
  217. % Change a picture v to {0,1} (`monochrome') form.
  218. % Apply this before the other monochrome operations.
  219.  
  220. def mono (suffix v) =
  221.   cull v keeping (1, infinity);
  222. enddef;
  223.  
  224. % bitwise and.
  225.  
  226. primarydef v picand w =
  227.  begingroup
  228.   picture u;  u=v+w;  cull u keeping (2,2);  u
  229.  endgroup
  230. enddef;
  231.  
  232. % inclusive or.
  233.  
  234. primarydef v picor w =
  235.  begingroup
  236.   picture u;  u=v+w;  cull u keeping (1,2);  u
  237.  endgroup
  238. enddef;
  239.  
  240. % exclusive or.
  241.  
  242. primarydef v picxor w =
  243.  begingroup
  244.   picture u;  u=v+w;  cull u keeping (1,1);  u
  245.  endgroup
  246. enddef;
  247.  
  248. % (nonsymmetric) difference.
  249.  
  250. primarydef v picsub w =
  251.  begingroup
  252.   picture u;  u=v-w;  cull u keeping (1,1);  u
  253.  endgroup
  254. enddef;
  255.  
  256. % contour to bitmap --- clipping and filling.
  257.  
  258. % (safely filled) interior of contour c,
  259. % where c is in graph coordinates;
  260. % adapted from _The MFbook_'s "safefill".
  261.  
  262. vardef interior expr c =
  263.   save vs; picture vs; vs=nullpicture;
  264.   interim turningcheck:=0;
  265.   addto vs contour (c.t_) withpen nullpen;
  266.   cull vs dropping (0,0);
  267.   vs
  268. enddef;
  269.  
  270. % derived bitmap operations.
  271.  
  272. % clip copy of picture vt to interior of cycle c,
  273. % where c is in graph coordinates,
  274. % return result (which is a `subpicture' of vt).
  275.  
  276. vardef clip (suffix vt) expr c =
  277.   mono (vt);
  278.   vt picand (interior c)
  279. enddef;
  280.  
  281. % fill region inside c in picture vt,
  282. % where c is in graph coordinates.
  283.  
  284. vardef picfill (suffix vt) expr c =
  285.   mono (vt);
  286.   vt := vt picor (interior c);
  287. enddef;
  288.  
  289. % unfill region inside c in picture vt,
  290. % where c is in graph coordinates.
  291.  
  292. vardef picunfill (suffix vt) expr c =
  293.   mono (vt);
  294.   vt := vt picsub (interior c);
  295. enddef;
  296.  
  297. % reverse video of vt inside c,
  298. % where c is in graph coordinates.
  299.  
  300. vardef picneg (suffix vt) expr c =
  301.   mono (vt);
  302.   (interior c) picsub vt
  303. enddef;
  304.  
  305. % rendering paths --- drawing and filling.
  306.  
  307. % draw path f in picture v using current pen,
  308. % where f is in graph coordinates.
  309.  
  310. def onepath (expr f) (suffix v) =
  311.  addto v doublepath (f.t_)
  312.    withpen currentpen;
  313. enddef;
  314.  
  315. % drawing a path d safely in picture vt,
  316. % where d is in graph coordinates.
  317. % (Courtesy of Uwe Bonnes.)
  318.  
  319. newinternal minpen;
  320. interim minpen := 0.01pt;
  321.  
  322. vardef picdraw (suffix vt) expr d =
  323.  save vs; picture vs; vs=nullpicture;
  324.  interim turningcheck:=0;
  325.  if penwd > minpen:
  326.    onepath (d, vs);
  327.    mono (vs);
  328.  fi;
  329.  addto vt also vs;
  330. enddef;
  331.  
  332. % drawing a path d safely,
  333. % where d is in graph coordinates.
  334.  
  335. def safedraw expr d =
  336.   picdraw (active_plane) d;
  337. enddef;
  338.  
  339. % filling a region, the interior of c, safely,
  340. % where c is in graph coordinates.
  341.  
  342. def safefill expr c =
  343.   if cycle c:
  344.     picfill (active_plane) c
  345.   else:
  346.     safedraw c
  347.   fi;
  348. enddef;
  349.  
  350. % erasing a region (the interior of c) safely,
  351. % where c is in graph coordinates.
  352. % - really works this time -
  353.  
  354. def safeunfill expr c =
  355.   if cycle c:
  356.     picunfill (active_plane) c
  357.   else:
  358.     safedraw c
  359.   fi;
  360. enddef;
  361.  
  362. % draw path f and return path f,
  363. % where f is in graph coordinates.
  364.  
  365. vardef drawnpath expr f =
  366.  safedraw f;
  367.  f
  368. enddef;
  369.  
  370.  
  371. %% Shading and Hatching.
  372.  
  373. % shading routine
  374. % in pixel coordinates
  375.  
  376. % one square dot scaled to (actual) size in x and y directions:
  377. % size is in pixel coordinates.
  378.  
  379. def setsqdot (expr size) =
  380.   if not picture sqdot:
  381.     picture sqdot;
  382.   fi
  383.   sqdot := nullpicture;
  384.   addto sqdot
  385.     contour unitsquare
  386.       xscaled (hround (size))
  387.       yscaled (vround (size))
  388.     withpen nullpen;
  389. enddef;
  390.  
  391. setsqdot (0.5pt);
  392.   
  393. % draw a square dot at position p in picture v,
  394. % where p is in _graph_ coordinates.
  395.  
  396. def squaredot (expr p)(suffix v) =
  397.   addto v also (sqdot shifted (hroundpair (p.t_)));
  398. enddef;
  399.  
  400. % calculate bounding box points (ll, ur) for path f.
  401. % This is coordinate system independent, but it's usually
  402. % employed when f is in pixel coordinates.
  403.  
  404. def bbox (expr f)(suffix ll, ur) =
  405.   ur := ll := point 0 of f;
  406.   pair p[];
  407.   for i=0 upto length f:
  408.     p0 := point i of f;
  409.     p1 := precontrol i of f;
  410.     p2 := postcontrol i of f;
  411.     ll := minpair (ll, p0, p1, p2);
  412.     ur := maxpair (ur, p0, p1, p2);
  413.   endfor;
  414.   % FOR TESTING ONLY: DRAW THE BOUNDING BOX:
  415.   % safedraw (rect (ll, ur));
  416. enddef;
  417.  
  418. % shade interior of path f with dots spaced sp apart,
  419. % where f is in _graph_ coordinates,
  420. % and sp is in _pixel_ coordinates.
  421.  
  422. def shadepath(expr sp) expr f =
  423.  if not cycle f:
  424.   safedraw f;
  425.  elseif sp<=0:
  426.   fill f;
  427.  else:
  428.   picture v;
  429.   pair p[], ll, ur, mn;
  430.   bbox (f, ll, ur);
  431.   ll:=sp*(ceilingpair(ll/sp));
  432.   mn:=floorpair((ur-ll)/sp);
  433.   m:=xpart mn;
  434.   n:=ypart mn;
  435.   twosp:=2*sp;
  436.   v:=nullpicture;
  437.   p2:=ll;
  438.   for i=0 upto m:
  439.    p3:=p2 if odd i: +(0,sp) fi;
  440.    for j=0 upto n:
  441.     if (not odd (i+j)):
  442.      % squaredot (p3, v);
  443.      onepath (p3, v);
  444.      p3:=p3+(0,twosp);
  445.     fi;
  446.    endfor;
  447.    p2:=p2+(sp,0);
  448.   endfor;
  449.   addto active_plane
  450.     also clip(v) f;
  451.  fi;
  452. enddef;
  453.  
  454. % hatch interior of path f with lines spaced sp apart,
  455. % where f is in _graph_ coordinates,
  456. % and sp is in _pixel_ coordinates.
  457.  
  458. def hatchpath(expr sp) expr f =
  459.  if not cycle f:
  460.   safedraw f;
  461.  elseif sp<=0:
  462.   fill f;
  463.  else:
  464.   picture v;
  465.   pair p[], ll, ur, mn;
  466.   bbox (f, ll, ur);
  467.   p0:=ll - (0,xpart ur - xpart ll);
  468.   p0:=sp*floorpair(p0/sp);
  469.   m:=ceiling((ypart ur + xpart ur
  470.     - ypart ll - xpart ll)/(sqrt(2)*sp));
  471.   v:=nullpicture;
  472.   p1:=(0,sqrt(2)*sp);
  473.   p2:=p0;
  474.   p3:=p2+(xpart ur - xpart ll + sp,
  475.     xpart ur - xpart ll + sp);
  476.   for i=0 upto m:
  477.    onepath (p2--p3, v);
  478.    p2:=p2+p1; p3:=p3+p1;
  479.   endfor;
  480.   addto active_plane
  481.     also clip(v) f;
  482.   p0:=ll + (xpart ur - xpart ll,
  483.     xpart ll - xpart ur);
  484.   p0:=sp*floorpair(p0/sp);
  485.   v:=nullpicture;
  486.   p2:=ll;
  487.   p2:=sp*floorpair(p2/sp);
  488.   p3:=p2+(xpart ur - xpart ll + sp,
  489.     xpart ll - xpart ur - sp);
  490.   for i=0 upto m:
  491.    onepath (p2--p3, v);
  492.    p2:=p2+p1; p3:=p3+p1;
  493.   endfor;
  494.   addto active_plane
  495.     also clip(v) f;
  496.  fi;
  497. enddef;
  498.  
  499.  
  500. %% Dots and Dashes.
  501.  
  502. % return estimated length of segment k of path f,
  503. % using an n-piece polyline approximation.
  504.  
  505. vardef seglength(expr f, k, n)=
  506.  pair p[], q[];
  507.  p1:=0.5[point k-1 of f,postcontrol k-1 of f];
  508.  p2:=0.5[postcontrol k-1 of f,precontrol k of f];
  509.  p3:=0.5[precontrol k of f,point k of f];
  510.  for i=2 upto n:
  511.   q1:=0.5[point k-1 of f,p1];
  512.   for j=2 upto i+1:
  513.    q[j]:=0.5[p[j-1],p[j]];
  514.    p[j-1]:=q[j-1];
  515.   endfor;
  516.   p[i+2]:=0.5[p[i+1],point k of f];
  517.   p[i+1]:=q[i+1];
  518.  endfor;
  519.  l:=length(p1-(point k-1 of f));
  520.  for i=1 upto n+1:
  521.   l:=l+length(p[i+1]-p[i]);
  522.  endfor;
  523.  l:=l+length((point k of f)-p[n+2]);
  524.  l
  525. enddef;
  526.  
  527. % draw dotted/dashed curve along path f,
  528. % with dash length dlen and dash space slen;
  529. % return f.
  530.  
  531. vardef dotted(expr dlen, slen) expr f =
  532.  pair p[]; numeric a[];
  533.  len:=0;
  534.  for i=1 upto length f:
  535.   a[i]:=seglength(f, i, 10);
  536.   len:=len+a[i];
  537.  endfor;
  538.  if len<dlen: safedraw f;
  539.  else:
  540.   n:=floor(len/(dlen+slen));
  541.   delta:=(len/n)-(dlen+slen);
  542.   dist:=0.5*dlen;
  543.   for i=0 upto (length f)-1:
  544.    len:=a[i+1];
  545.    forever:
  546.     exitif dist>len;
  547.     safedraw subpath
  548.       (i if dist>dlen: +(dist-dlen)/len fi, i+(dist/len))
  549.       of f;
  550.     dist:=dist+slen+delta+dlen;
  551.    endfor;
  552.    dist:=dist-len;
  553.    if dist<dlen:
  554.     safedraw subpath (i+(len+dist-dlen)/len, i+1) of f;
  555.    fi;
  556.   endfor;
  557.  fi;
  558.  f
  559. enddef;
  560.  
  561.  
  562. %% Points.
  563.  
  564. def pointd(expr ptwd,filled)(text t) =
  565.  pair p.px;
  566.  for a=t:
  567.   p.px:=a transformed ztr;
  568.   if filled: safefill else: safeunfill drawnpath fi
  569.     fullcircle scaled ptwd shifted p.px;
  570.  endfor
  571. enddef;
  572.  
  573.  
  574. %% Arrows.
  575. % arrowheads are in pixel coordinates
  576.  
  577. boolean hfilled;
  578. newinternal hdwdr, hdten;
  579. interim hdwdr:=1;
  580. interim hdten:=1;
  581.  
  582. def head(expr front, back, width, t, filled) =
  583.  pair p[], side;
  584.  side := (width/2) *
  585.    ((front-back) rotated 90);
  586.  p1 := back + side;
  587.  p2 := back - side;
  588.  if filled: safefill drawnpath else: safedraw fi
  589.    p1..tension t..{front-back}(front){back-front}..tension t..p2
  590.    if filled: --cycle; fi;
  591. enddef;
  592.  
  593. vardef headpath(expr hlen, hrot, hback) expr f=
  594.  pair p[];
  595.  p2:=point length f of f;
  596.  p1:=direction length f of f;
  597.  if p1<>(0,0):
  598.   p3:=(unitvector(p1) rotated hrot);
  599.   head(p2-(hback*p3),p2-((hlen+hback)*p3),
  600.     hdwdr,hdten,hfilled);
  601.  fi;
  602.  f
  603. enddef;
  604.  
  605. def arrow(expr tl,hd,hlen) =
  606.  path f.px;
  607.  f.px:= (tl..hd) transformed ztr;
  608.  safedraw f.px;
  609.  trash:=headpath(hlen,0,0pt) f.px;
  610. enddef;
  611.  
  612.  
  613. %% Axes and Axis Tic Marks.
  614.  
  615. def axes(expr hlen) =
  616.  arrow((0,yneg),(0,ypos),hlen);
  617.  arrow((xneg,0),(xpos,0),hlen);
  618. enddef;
  619.  
  620. def xmarks(expr len)(text t) =
  621.  for a=t:
  622.   safedraw (xconv(a),yconv(0)-(len/2))..
  623.     (xconv(a),yconv(0)+(len/2));
  624.  endfor;
  625. enddef;
  626.  
  627. def ymarks(expr len)(text t) =
  628.  for a=t:
  629.   safedraw (xconv(0)-(len/2),yconv(a))..
  630.     (xconv(0)+(len/2),yconv(a));
  631.  endfor;
  632. enddef;
  633.  
  634.  
  635. %% Path Construction.
  636.  
  637. % polyline construction
  638.  
  639. def mkpoly (expr cyclic)
  640.   (suffix pts) =
  641.  for i=1 upto pts-1:
  642.   pts[i]--
  643.  endfor
  644.  pts[pts]
  645.  if cyclic:
  646.   --cycle
  647.  fi
  648. enddef;
  649.  
  650. % smooth path construction
  651.  
  652. def mksmooth (expr cyclic)
  653.   (suffix pts) =
  654.  if cyclic:
  655.   pts[1]{pts[2]-pts[pts]}
  656.  else:
  657.   pts[1]
  658.  fi
  659.  for i=2 upto pts-1:
  660.   ..pts[i]{pts[i+1]-pts[i-1]}
  661.  endfor
  662.  if cyclic:
  663.   ..pts[pts]{pts[1]-pts[pts-1]}..cycle
  664.  else:
  665.   ..pts[pts]
  666.  fi
  667. enddef;
  668.  
  669. % general path construction
  670.  
  671. def mkpath(expr smooth, cyclic)
  672.   (suffix pts) =
  673.  if smooth:
  674.   mksmooth(cyclic,pts)
  675.  else:
  676.   mkpoly(cyclic,pts)
  677.  fi
  678. enddef;
  679.  
  680.  
  681. %% Upright Rectangles.
  682.  
  683. def mkrect(expr ll,ur) =
  684.  ll--(xpart ll,ypart ur)--
  685.    ur--(xpart ur,ypart ll)--cycle
  686. enddef;
  687.  
  688. def rect(expr ll,ur) =
  689.  (mkrect(ll,ur))
  690. enddef;
  691.  
  692.  
  693. %% Curves.
  694.  
  695. vardef mkcurve(expr smooth,cyclic)
  696.   (text t) =
  697.  textpairs(p,t);
  698.  mkpath(smooth,cyclic,p)
  699. enddef;
  700.  
  701. def curve(expr smooth,cyclic)
  702.   text t =
  703.  mkcurve(smooth,cyclic,t)
  704. enddef;
  705.  
  706. % quadratic B-splines.
  707. % p[] == B-spline control points,
  708. % p in number.
  709.  
  710. vardef openqbs (text t) =
  711.   textpairs(p,t);
  712.   for i=1 upto p-2:
  713.     0.5[p[i],p[i+1]]
  714.     ..controls 1/6[p[i+1],p[i]]
  715.            and 1/6[p[i+1],p[i+2]]..
  716.   endfor
  717.   0.5[p[p-1],p[p]]
  718. enddef;
  719.  
  720. vardef closedqbs (text t) =
  721.   textpairs(p,t);
  722.   p[p+1]:=p1;
  723.   p[p+2]:=p2;
  724.   for i=1 upto p:
  725.     0.5[p[i],p[i+1]]
  726.     ..controls 1/6[p[i+1],p[i]]
  727.            and 1/6[p[i+1],p[i+2]]..
  728.   endfor
  729.   cycle
  730. enddef;
  731.  
  732. % cubic B-splines.
  733.  
  734. vardef mkopencbs (suffix b) =
  735.   for i=1 upto b-3:
  736.     (b[i]+4b[i+1]+b[i+2])/6
  737.     ..controls 1/3[b[i+1],b[i+2]]
  738.            and 2/3[b[i+1],b[i+2]]..
  739.   endfor
  740.   (b[b-2]+4b[b-1]+b[b])/6
  741. enddef;
  742.  
  743. vardef opencbs (text t) =
  744.   textpairs(b,t);
  745.   mkopencbs(b)
  746. enddef;
  747.  
  748. vardef mkclosedcbs (suffix b) =
  749.   b[b+1]:=b1;
  750.   b[b+2]:=b2;
  751.   for i=1 upto b:
  752.     (b[i]+4b[i+1]+b[i+2])/6
  753.     ..controls 1/3[b[i+1],b[i+2]]
  754.            and 2/3[b[i+1],b[i+2]]..
  755.   endfor
  756.   cycle
  757. enddef;
  758.  
  759. vardef closedcbs (text t) =
  760.   textpairs(b,t);
  761.   mkclosedcbs(b)
  762. enddef;
  763.  
  764.  
  765. %% Path Closure.
  766.  
  767. % path f closed by line segment.
  768.  
  769. vardef closedpath expr f =
  770.  f
  771.  if not cycle f:
  772.    --cycle
  773.  fi
  774. enddef;
  775.  
  776. % close path f in manner of mksmooth.
  777.  
  778. vardef sclosedpath expr f =
  779.  if cycle f:
  780.   f
  781.  else:
  782.   save n;
  783.   n := length f;
  784.   if n >= 2:
  785.    (point 0 of f){(point 1 of f)-(point infinity of f)}
  786.    ..(subpath (1,n-1) of f)
  787.    ..(point infinity of f){(point 0 of f)-(point (n-1) of f)}
  788.    ..cycle
  789.   elseif n = 1:
  790.    f--cycle
  791.   else: % single point
  792.    f..cycle
  793.   fi
  794.  fi
  795. enddef;
  796.  
  797. % path f closed by bezier.
  798.  
  799. vardef bclosedpath expr f =
  800.  f
  801.  if not cycle f:
  802.    ..cycle
  803.  fi
  804. enddef;
  805.  
  806. % conversion of Bezier segment key points, z,
  807. % to cubic B-spline control points, b.
  808.  
  809. def ztob (suffix z, b) =
  810.   pair b[];
  811.   b := 4;
  812.   b1 = 6z1-7z2+2z3;
  813.   b2 =     2z2- z3;
  814.   b3 =    - z2+2z3;
  815.   b4 =     2z2-7z3+6z4;
  816. enddef;
  817.  
  818. % closure of path f by a cubic B-spline.
  819.  
  820. vardef cbclosedpath expr f =
  821.  if cycle f:
  822.   f
  823.  else:
  824.   pair p[];
  825.   p1 := point 0 of f;
  826.   p2 := postcontrol 0 of f;
  827.   p3 := precontrol 1 of f;
  828.   p4 := point 1 of f;
  829.   ztob(p,a);
  830.   n := length f;
  831.   p1 := point (n-1) of f;
  832.   p2 := postcontrol (n-1) of f;
  833.   p3 := precontrol infinity of f;
  834.   p4 := point infinity of f;
  835.   ztob(p,b);
  836.   b1 := b3;
  837.   b2 := b4;
  838.   b3 := a1;
  839.   b4 := a2;
  840.   f..mkopencbs(b)..cycle
  841.  fi
  842. enddef;
  843.  
  844.  
  845. %% Circles and Ellipses.
  846.  
  847. vardef mkellipse(expr center,radx,rady,
  848.   angle) =
  849.  save t;
  850.  transform t;
  851.  t:=identity xscaled (2*radx)
  852.    yscaled (2*rady) rotated angle
  853.    shifted center;
  854.  fullcircle transformed t
  855. enddef;
  856.  
  857. def ellipse(expr center,radx,rady,
  858.   angle) =
  859.  (mkellipse(center,radx,rady,angle))
  860. enddef;
  861.  
  862. def circle(expr center,rad) =
  863.  ellipse(center,rad,rad,0)
  864. enddef;
  865.  
  866.  
  867. %% Circular Arcs.
  868.  
  869. vardef mkarc(expr center,from,sweep)=
  870.  pair p,q;
  871.  path f;
  872.  if sweep=0: f:=from
  873.  else:
  874.   n:=floor(abs(sweep)/45)+1;
  875.   if n<3: n:=3; fi;
  876.   theta:=sweep/(n-1);
  877.   f:=p:=from;
  878.   for i=2 upto n:
  879.    p:=p rotatedabout (center,theta);
  880.    q:=p-center; q:=q rotated 90;
  881.    if theta<0: q:=-q; fi;
  882.    f:=f..p{unitvector q};
  883.   endfor;
  884.  fi;
  885.  f
  886. enddef;
  887.  
  888. vardef arccenter(expr from,to,sweep)=
  889.  pair midpt;
  890.  if from=to:
  891.    from
  892.  else:
  893.    midpt:=(0.5)[from,to];
  894.    if (sweep mod 360)=0:
  895.      midpt
  896.    else:
  897.      disp:=cosd(sweep/2)/sind(sweep/2);
  898.      midpt+(disp*((to-from) rotated 90)/2)
  899.    fi
  900.  fi
  901. enddef;
  902.  
  903. % point-point-sweep form of arc
  904.  
  905. vardef arcpps(expr from,to,sweep) =
  906.  pair center;
  907.  center:=arccenter(from,to,sweep);
  908.  (mkarc(center, from, sweep))
  909. enddef;
  910.  
  911. % modified polar --
  912.  
  913. % center, angle, angle, radius
  914.  
  915. vardef arcplr(expr center,
  916.   frtheta,totheta,rad) =
  917.  pair from;
  918.  from:=center+rad*(dir frtheta);
  919.  (mkarc(center,from,
  920.    totheta-frtheta))
  921. enddef;
  922.  
  923. % center-point-sweep form
  924.  
  925. def arccps(expr center, from, theta)=
  926.  (mkarc(center, from, theta))
  927. enddef;
  928.  
  929. % point-point-point form
  930.  
  931. vardef arcppp(expr first, second, third)=
  932.  sweep:=2*(angle(third-second)-angle(second-first));
  933.  sweep:=sweep mod 720;
  934.  if sweep > 360: sweep:=sweep-720; fi
  935.  critical:=5;
  936.  if abs(sweep) <= critical:  % center may blow out
  937.   save p;
  938.   pair p[];
  939.   p:=3;
  940.   p1:=first;
  941.   p2:=second;
  942.   p3:=third;
  943.   mkpath(true,false,p)
  944.  else:
  945.   pair m[], d[], center;
  946.   m1:=(0.5)[first,second]; d1:=(second-first) rotated 90;
  947.   m2:=(0.5)[second,third]; d2:=(third-second) rotated 90;
  948.   center = m1+whatever*d1 = m2+whatever*d2;
  949.   mkarc(center,first,sweep)
  950.  fi
  951. enddef;
  952.  
  953.  
  954. %% Polar Coordinates.
  955.  
  956. % (r, theta) -> (x, y).
  957.  
  958. def plrconv(expr a)=
  959.  ((ypart a)*(dir xpart a))
  960. enddef;
  961.  
  962. def plrpointd(expr ptwd,filled)(text t) =
  963.  pair p.px;
  964.  for a=t:
  965.   p.px:=plrconv(a) transformed ztr;
  966.   if filled: safefill else: safeunfill drawnpath fi
  967.     fullcircle scaled ptwd shifted p.px;
  968.  endfor
  969. enddef;
  970.  
  971. vardef mkplrcurve(expr smooth,cyclic)
  972.   (text t)=
  973.  save p;
  974.  pair p[];
  975.  p:=0;
  976.  for a=t:
  977.   p[incr p]:=plrconv(a);
  978.  endfor;
  979.  mkpath(smooth,cyclic,p)
  980. enddef;
  981.  
  982. def plrcurve(expr smooth,cyclic)
  983.   text t =
  984.  mkplrcurve(smooth,cyclic,t)
  985. enddef;
  986.  
  987. % other figures
  988.  
  989. vardef turtle(text t)=
  990.  save p;
  991.  pair p[];
  992.  p:=0;
  993.  th:=0;
  994.  for a=t:
  995.   if p=0:
  996.    p[incr p]:=a;
  997.   else:
  998.    th:=th+(xpart a);
  999.    p[incr p]:=((ypart a)*(dir th));
  1000.    p[p]:=p[p]+p[p-1]-((0,0) );
  1001.   fi;
  1002.  endfor;
  1003.  (mkpath(false,false,p))
  1004. enddef;
  1005.  
  1006.  
  1007. %% Wedges.
  1008.  
  1009. vardef wedge(expr center,
  1010.   frtheta,totheta,rad) =
  1011.  pair from;
  1012.  from:=center+rad*(dir frtheta);
  1013.  (center--
  1014.    mkarc(center,from,totheta-frtheta)
  1015.    --cycle)
  1016. enddef;
  1017.  
  1018.  
  1019. %% Functions.
  1020.  
  1021. vardef mkfcn(expr smooth,bmin,bmax,bst)
  1022.   (suffix bv)(text fcnpr)=
  1023.  save p;
  1024.  pair p[];
  1025.  p:=0;
  1026.  for bv=bmin step bst until bmax+(bst/2):
  1027.   p[incr p]:=fcnpr;
  1028.  endfor;
  1029.  mkpath(smooth,false,p)
  1030. enddef;
  1031.  
  1032. def function(expr smooth,xmin,xmax,st)
  1033.   (text fx) =
  1034.  (mkfcn(smooth,xmin,xmax,st,
  1035.    x,(x,fx)))
  1036. enddef;
  1037.  
  1038. def parafcn(expr smooth,tmin,tmax,st)
  1039.   (text ft) =
  1040.  (mkfcn(smooth,tmin,tmax,st,
  1041.    t,ft))
  1042. enddef;
  1043.  
  1044. def plrfcn(expr smooth,tmin,tmax,st)
  1045.   (text ft) =
  1046.  (mkfcn(smooth,tmin,tmax,st,
  1047.    t,((ft)*cosd(t),(ft)*sind(t))))
  1048. enddef;
  1049.  
  1050. %%%
  1051. %%%  end  graphbase.mf
  1052. %%%
  1053.